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ABSTRACT 

We use high resolution N-body simulations to investigate the Lagrangian bias of cold 
dark matter haloes within the A cold dark matter cosmology. Our analysis focuses on 
'protohaloes', which we identify in the simulation initial conditions with the subsets of 
particles belonging to individual redshift-zero haloes. We then calculate the number- 
density and velocity-divergence fields of protohaloes and estimate their auto spectral 
densities. We also measure the corresponding cross spectral densities with the linear 
matter distribution. We use our results to test a Lagrangian-bias model presented 
by Desjacques and Sheth which is based on the assumption that haloes form out of 
local density maxima of a specific height. Our comparison validates the predicted 
functional form for the scale-dependence of the bias for both the density and velocity 
fields. We also show that the bias coefficients arc accurately predicted for the velocity 
divergence. On the contrary, the theoretical values for the density bias parameters do 
not accurately match the numerical results as a function of halo mass. This is likely 
due to the simplistic assumptions that relate virialized haloes to density peaks of a 
given height in the model. We also detect appreciable stochasticity for the Lagrangian 
density bias, even on very large scales. These are not included in the model at leading 
order but correspond to higher order corrections. 

Key words: mcthods:analytical - numerical - galaxies: haloes - cosmology: theory 
- large-scale structure of Universe. 
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1 INTRODUCTION 

Galaxy redshift surveys are powerful probes of cosmology. 
The main observables able to constrain cosmological param- 
eters are the overall shape of the galaxy power spectrum 
at wavenumbers k < 0.1 Mpc -1 and the baryonic acoustic 
oscillations within it. These are treated as proxies for the 
matter power spectrum for which we can make robust the- 
oretical predictions. Galaxies, however, are biased tracers 
of the cosmic mass distribution and many features appear- 
ing in their power spectrum depend on how a specific ob- 
servational sample was selected. To reconstruct the matter 
power spectrum we thus need an accurate bias model whose 
free coefficients should be used as nuisance parameters and 
marginalized over. In the era of precision cosmology, where 
measurements of the matter power spectrum with per cent 
accuracy are required, this task is particularly demanding. 

Bias models can be divided into two broad classes. Eu- 
lerian biasing schemes relate the galaxy density contrast, 
<S g (x,t), to the matter density distribution, S, evaluated 
at the same time t (but not necessarily at the same spa- 
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tial location). After smoothing the fields on large scales, so 
that \S\ is typically muc h smaller than unity, one can write 
(|Frv fc Gaztanagalll993h 



<5 g (x) = B + / d 3 a:iBi(x-xi)5(x 1 ) + 



+ 
+ 



(1) 



/ 3 3 

- / d xi d x 2 B 2 (x — xi, x — x 2 ) <S(xi) <5(x 2 ) + 



where all fields are evaluated at the same time t and the 
details of the bias model are specified by the kernel func- 
tions, Bi. If one further assumes that biasing is local (i.e. 
that all kernels can be written as products of Dirac delta 
distributions), this reduces to 



5 (x) = 6 + 6i5(x) + ^ 2 (x) + 



(2) 



where now the bias coefficients bi are real numbers. 

In Lagrangian bias models, on the other hand, one con- 
siders the regions in the initial conditions that will collapse 
to form galaxies (or their hosting dark-matter haloes) at 
time t and writes their density contrast, <5g(q), in terms 
of the linear density contrast, 5o(q)- Large-scale expansions 
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analogous to eqs. ([T]) and ([2]) can also be written in this case. 
As a second step, one must determine the final position, 
x(q, t), of a fluid element initially located at q, and com- 
pute 5 g (x(q, t),t) out of 5g(q). This Lagrangian-to-Eulerian 
mapping (LEM) accounts for gravitationally induced mo- 
tions that determine the final position of the objects. 

Local Eulerian and Lagrangian bias schemes are not 
equivalent; t hey generate a diffe rent shape for the galaxy 
bispectrum (|Catelan et alj l2000h and are not compatible 
withi n the fr amework of perturba tion theory (jMatsubaral 
l201ll ). In fact. ICatelan et all (|l998T ) have shown that a local 
Lagrangian biasing scheme generates a non- linear, non-local 
and stochastic bias in Eulerian space. On the other hand, 
non-local Eulerian and Lagrangian schemes are equivalent 
and can be seen as different mathematical re presentations 
of the same physical process (jMatsubaral |20 1 lT ) . 

Due to its simplicity, the local Eulerian model is by 
far the most widely used in practical applications, such 
as perturbative calculations. However, it is purely phe- 
nomenological and does not have a strong theoretical mo- 
tivation. Detailed comparison with num erical simulations 
has a l so evidenced its limited validity (e.g. iRoth fe Porcianil 
l201ll . IPollack. Smith fc Porcianil [201lh . Physical models of 
bias are generally given in the Lagrangian framework as 
conditions on galaxy (halo) formation are more easily 
imposed onto the linear density field using some model 
for th e collapse of density perturbations. iMo fc White! 
19961. hereafter MW) use d a Press-Schechter-like argument 



Press fc Schechter1ll974D to compute the bias coefficients 



of a local Lagrangian scheme as a function of halo mass. 
The same authors also showed how these parameters can be 
combined to calculate the bias coefficients of a local Eule- 
rian scheme assuming that large-scale density perturbations 
follow the spherical collapse model. T he effect of non-linear 
shear on the LEM was discussed by ICatelan et al. Jl998T ) 
using the Zel'dovich approximation ljZerdovichlll97C ) . For 
halo masses M > M*(z) (where M,(z) is the character- 
istic mass for collapse at redshift z) the MW formula for 
fei is in goo d agreement with t he pr edictions of N-body 
simulations (|Mo. Jing fc White! Il997f ). At lower masses, 
however, the agreement r a pidly deter iorates jjind Il998f ) . 
IPorciani. Catelan fc Lacevl (|l999h and Ijind l| 1999h showed 
that the discrepancy between the N-body simulations and 
the analytical predictions is already present in Lagrangian 
space and should thus be attributed to the limitations of the 
Press-Schechter formalism rather than to the approximated 
treatment of the LEM. 

The Lagrangian bias em erging from the ex- 
tended Press-Sch echter model dBond et al.l Il99l[ ) was 
_ Porciani et al.l ( 1998h . rediscussed 



first 



derived by „ 

IScannapieco fc Barkanal J2002h an d teste d against simu- 
lations by IScannapieco fc Thackerl (|2005h . This approach 
follows correlated trajectories of Sq at different Lagrangian 
locations as a function of the smoothing scale and looks 
for correlations in the first-crossing scales of a density 
threshold. 

According to the peak-backgroun d-split argument 
(jBardeen et al.lll986l : ICole fc Kaiserlll989l ). long- wavelength 
density fluctuations modulate halo formation by modifying 
the collapse time of localized short-wavelength perturba- 
tions. This makes it possible to generalise the calculation 
of the Lagrangian MW bias coefficients to any model for 



the halo mass function dSheth fc Tormenlll999l ,l Tinker et al.l 
2010. Giannantonio fc Porcianil |2010|) and to improve the 
agreement with N-body simulations Jseliak fc Warren! [2004I . 



I Tinker et al 



Kaiser 



20051 . iGao et al.ll2005l . iPillepich et al.ll2010T ) 



1984T ) explained the strong clustering of Abell 



clusters by assuming that they originate from the regions 
above a density threshold in the (suitably smoothed) lin- 
ear density field. Following this line of reasoning, it is com- 
mon to assume that dark-matter haloes form out of lin- 
ear density peaks, as an alternative to the Press-Schechter 
approach. Tests against N-body simulations have shown 
that this is a well justified assumptio n, especi ally for mas- 
sive haloes ([Ludlow fc Porcianil l201ll . see also iFrenk et al.l 
Il988f l. The statistical properties of the local maxima in a 
Gaussian random fi eld have been extensively studied by 
iBardeen et al.l (Tl986h (see a lsolPeacock fc H eavens 1985| and 
iHoffman fc Shahamlfl985h . IMo. Jing fc White! Jl997lnntro- 



duced peaks theory in the MW formalism, while iMatsubaral 
(jl999f l evaluated the level of stochasticity in the Lagrangian 

clustering of d ensity extrema. 

Recently, iDesiacauesl l|2008f ) and iDesiacaues fc Shethl 

|2010l . hereafter DS) showed that the correlation of den- 
sity peaks in real and redshift space can be interpreted in 
terms of a simple Lagrangian biasing scheme. Due to the 
peak constraint, the effective Lagrangian peak density at 
a given point not only depends on the local value of the 
mass density but also on its Laplacian. The first-order peak 
bias depends on the mass and height of the peaks, and on 
the matter po wer spectrum. For high peaks, this reduces to 
the results bv IMatsubaral ll 19991 ). Moreover, although peaks 
move with the matter at their positions, DS infer the exis- 
tence of a statistical velocity bias due to the fact that local 
maxima can only exist at special locations. This also leads 
to a bias between the linear velocity spectra of peaks and 
matter which is predictable in quantitative terms. 

In spite of the fact that the DS model provides the ini- 
tial conditions for sophisticated models of the (Eulerian) 
halo distribution where the L EM is based ei ther on re- 
summed perturbation the ory jElia et all 2011) o r on the 
Zel'dovich approximation ijDesiacaues et al .1 120101 ) . its pre- 
dictions for the Lagrangian clustering and velocities of the 
regions that will form collapsed structures have never been 
thoroughly tested against numerical simulations. This pa- 
per provides such a test, which is necessary if we are to use 
advanced bias models to extract useful information on the 
cosmological parameters through a comparison with obser- 
vations. 

The structure of the paper is as follows. In Section [2] 
we review the bias model first presented in DS. The details 
of our N-body simulations and analysis techniques are out- 
lined in Section[3l and a comparison between the model and 
numerical results is presented in Section [4] Finally, we sum- 
marize our main results in Section [5] 



2 THE DS MODEL 

In this section we summarize, for completeness, the peaks 
model described in Section 2 of DS. In order to do so, we 
first introduce and define some quantities relevant for peak 
statistics. 
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The spectral moments of the matter power spectrum 
are defined as 



and 



27T 2 



_2,r> _\ _ 1 / J7.i,2(>>+1) r>/i. -Mirti. d \2 



Akk^ n+L> P(k,z)W(k,R s 



(3) 



where P(k, z) is the linear matter power spectrum at red- 
shift z, and W(k, R s ) is a smoothing kernel of characteristic 
length R s . In terms of the moments we define the spectral 
parameters: 

7n = a "I ■ (4) 

Un-lCTn+l 

There are two common choices for W{k,R s ): the Gaussian 
filter 

W G (k,R B ) , (5) 



and the top-hat filter 

xir n n \ 3[sin(fcP s ) - kR s cos(fcP B )] 
W T u(k,R s ) = -L-i ^3 i — a . 

The mass contained within the comoving length R s in these 



(6) 



M G (R S ) = (27r) 3/2 p7? a 3 , and M TH (Ps 



(7) 



where p is the mean matter density of the Universe. 

We will often characterize peaks in terms of their di 
mensionless peak height, v, defined 

S,:, 



v(Rs, z c ) = 



(8) 



ao(R 8 , 2c) ' 

Here 5 P is the smoothed overdensity at the peak location lin- 
early extrapolated to z = 0, and <jq(R s , z c ) is the linear rms 
mass fluctuation in spheres of radius Rs. DS linked density 
maxima of height v(Rs, z) to dark matter haloes of mass M s 
collapsing at redshift z c , assuming that <5 P coincides with the 
thre shold for collap s e, S c - 

iBardeen et all (|l986h and DS computed the cross- 
correlation between peaks and the underlying density field 
which also corresponds to t he average density profile around 
density maxima. Similarly, IPesiacquesI (|2008r i evaluated the 
leading order expressions (on large spatial separations) for 
the peak auto-correlation function and the line-of-sight 
m ean streaming for p airs of discrete local maxima of height 
v. IPesiacquesI (|2008l ) and DS showed that the full expres- 
sion for the cross-correlation and the large-scale asymptotic 
of the auto-correlation function are consistent with - and 
thus can be thought of as arising from - an effective bias 
relation. In the DS model, the number density and the ve- 
locity of peaks, <5n p k and v p k, are related to the dark matter 
density contrast and velocity fields, linearly extrapolated to 
z = 0, via: 

<5n pk (x|z c ) = £v<5s(x) - i> c V 2 <5s(x), (9) 



and 



<T0 2 , 



Vpk(x|;z c ) = Vs(x) j V(5s(x) 



(10) 



The subscript "S" indicates that the fields are smoothed on 
the scale R s , and the bias parameters, 6„ and , are given 
by 



l 

00 



7i u 



Of 



(11) 



<?2 



U — J\V 



(12) 



Here u is the me an curvature of the peaks, which can be 
approximated by (|Bardeen et al.|[l986l ) 



3(1 - 7l 2 ) + (1.216 - exp 



u = 7i^+- 



71 / Tig ' 
2 V 2 j 



3(1 - 7l 2 ) + 0.45 



+ (T) 



•(13) 



Note that b v coincid es with the peak bias factor found by 
IBardeen et al.l (|l986h after neglecting the derivatives of the 
density correlation function. 

Since, by definition, the gradient of the density field 
vanishes at peak locations, eq. (|10p suggests that the peak 
and dark matter velocity fields must be coincident there. By 
construction, peaks move with the dark matter flow, yet the 
spatial bias induces a statistical velocity bias. We will con- 
sider the scaled velocity divergence, #(x) = Vxv(x)/(aF/), 
rather than the velocity field. Here a is the scale factor, H 
the Hubble parameter, / = din D/d\na, with D the linear 
growth factor. In these units, both #(x) and <5(x) are dimen- 
sionless quantities. With these changes, eqs. Q and (|10|l 
can be rewritten in Fourier space as 



<5n pk (k) = (K + h c k 2 ) <5(k) W{k, R s 



and 



Mk) 



b a k 2 ) 0(k) W(k, Rs) = b vcl (k)6{k) . 



where we have defined 



b -d 



(14) 



(15) 



(16) 



In the limit of high peaks [y 3> 1) it can be shown that 
the bias parameters obey the following asymptotic relations: 
b v — > v/ctq and b( 0. This implies that the highest 
peaks are linearly biased tracers of the underlying matter 
field. This is con sistent with the predicti ons of the peak- 
background split (|Mo. Jing fc Whit j[l997l) and DS showed 
that, indeed, b v is the appropriate generalization of the con- 
stant, large-scale bias for low v. Unlike the density bias fac- 
tors, b a does not depend on v. 

In order to test this model against N-body simulations, 
we will make use of the cross-spectra between the peak 
and dark matter densities and velocities (denoted as P mp 
and Vmp, respectively) and of the corresponding peak auto- 
spectra (P p and P p ). From eqs. ()14|) and (|15[) we obtain 

Pnp(fc) = {b„ + b ( k 2 ) P{k) W{k,R s ), 



Pp(fc) ~ (K +b ( k 2 ) 2 P(k)W 2 {k,R s ) 

P mp (fc) ~ (l-Kk 2 )r(k)W{k,R B ), 
Pp(fc) ~ (l-b„k 2 ) 2 P(k)W 2 (k,R s ) 



(17) 



(18) 



where P(k) and V{k) are the matter density and veloc- 
ity divergence auto-spectra, respectively. We remind the 
reader that the expression for P P (k) is only valid to first 
order in P(k) as k — > 0, and that higher order correc- 
tions should be includ ed to improve its accuracy (see e.g. 
iDesiacques et al]|2010l ). On the contrary, the expression for 
the cross-spectrum P mp is exact, as shown in IBardeen et al.l 
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(| 19861 ) and in the Appendix A of DS. Eg dT71) h as the same 
functional form as eq. (57) in IMatsubaral (|l999h who stud- 
ied the clustering of de nsity extrema. The DS coefficients b v 
and match those in IMatsubaral 1 19991 ) only in the limit 
c>l, for which nearly all extrema are density maxima. 



3 NUMERICAL ISSUES 

In this section we provide a brief description of the main 
numerical issues relevant for this work. This includes a brief 
description of our numerical simulations in Section f3. II our 
main analysis techniques in Section T3. 21 and a characteriza- 
tion of halo collapse barriers in Section 13.31 



3.1 N-body simulations 

Our analysis focuses on two high-resolution N-body simu- 
lations of structure formation in the standard LCDM cos- 
mology. The cosmological parameters for our runs were cho- 
sen to be consistent w ith the fifth- year WMAP data release 
(|Komatsu et al.l 120091 '). These are h = 0.701, <r 8 = 0.817, 
n s = 0.96, n m = 0.279, ft b = 0.0462 and fi A = 1 - fim = 
0.721. Each simulation wa s run with a l ean version of the 
Tree-PM code Gadget-2 (|Springelll20"05T ), and followed the 
dark matter using 1024 3 collisionless particles. One simula- 
tion had a box side-length of Lbox = 1200 ft -1 Mpc and a 
particle mass of m par t = 1.246 ■ 10 11 h -1 M@; the other used 
L box = 150 h' 1 Mpc and had m part = 2.433 ■ lO 8 h -1 M . 
The initial redshifts of the simulations were Z[ n — 50 and 
Zi n = 70 for the larger and smaller box, respectively. Using 
these simulations we are able to probe a wide range of halo 
masses, spanning 8- lO lo /i _1 M < M h < h^My. These 
simulati ons were first studied in | Pillepich et al.l (|2010h , and 



later bv lLudlow fc Porcianil (|2011 ). and we refer the reader 



to those papers for further details. 

Haloes were i dentified at z = using a friends-of-friends 
(FOF) algorithm |Davis et al.ll 19861 ) with a linking length of 
0.2 times the mean interparticle distance. Protohaloes were 
identified by tracing backward to the initial redshift of all 
subsets of the particles belonging to z = FOF haloes. We 
use the centre of mass of each protohalo as a proxy for its 
spatial location; the mass-weighted linear velocity provides 
an estimate of the protohalo's motion. As a test of the sen- 
sitivity of our results to the adopted halo finder, we also 
generated a spherical overdensity (SO) halo catalogue with 
an overdensity threshold of 200 times the critical density, 
p c . For a fixed halo mass, the two halo-finders produce re- 
sults consistent within 10 per cent (in terms of all the bias 
coefficients) and so, in what follows, we will focus on re- 
sults obtained for the FOF haloes, and only consider those 
containing at least 100 particles. 

Haloes in each simulation are split into four separate 
mass bins in order to preserve their peculiar clustering prop- 
erties. These bins are referred to as IS to 4S for the small 
box, and 1L to 4L for the large one. To asses the impact 
of shot noise in the analysis of our small-box simulation we 
consider an additional mass bin, labeled binOS, which in- 
cludes all haloes with N > 100. The mass ranges and total 
number of haloes in each bin are given in Table [T] We note 
that bins with label 'S' refer to masses M < M* for which 



dark-matter haloes are not exp ected to be in one-to-one c or- 
respondence with linear peaks |Ludlow fc Porcianill201 J ). 



3.2 Analysis 

We construct protohalo density and momentum fields us- 
ing cloud-in-cell grid assignment on a 512 3 mesh. Velocity 
fields are obtained by taking the ratio of the momen tum 
and density fields, as described in IScoccimarrol (|2004l ). In 
the case of haloes, these distributions are smoothed to pre- 
clude the existence of empty cells; the smoothing scales used 
are Rf = 7h~ Mpc for the large box and R{ = 1.8ft Mpc 
for the small one. These values are chosen in order to mask 
the effects of the grid, but we have explicitly verified that 
our results are not significantly affected by them. All power 
spectra have been computed using a fast Fourier transform 
technique. 

The discreteness of dark-matter particles and haloes 
gives rise to a shot-noise component in the spectra. For the 
density fields, the estimated power spectrum, P, includes a 
shot-noise term which is inversely proportional to the num- 
ber density of objects n (assuming Poisson sampling): 



P — ftrue + 



(19) 



Shot noise is therefore negligible for the matter spectra but 
may be significant for that of the haloes. The issue is more 
severe in Lagrangian space, because fluctuations in the ini- 
tial conditions are small. We will consider two alternative 
estimates of the protohalo bias; one is determined from the 
shot-noise corrected auto-spectrum, 



b{k) 



Ph(fc) 
P(k) 



(20) 



and the other from the cross-spectrum, 

Here the subscript "h" indicates the halo fields, and "m" 
the matter field (the analogous fields for peaks are indicated 
with the subscript "p".) The relation between the two is 



6 eff (fe) = b(k) • r(k) , 

where r is the linear correlation coefficient, defined as 



r(k) 



Pmh(fc) 

y/mPjk) 



(22) 



(23) 



These relations tell us that the two definitions of the bias are 
equivalent only if the bias is purely deterministic in Fourier 
space, i.e. r = 1. At leading order, in the model presented 
by DS, fo c ff (k) — b v + b^ k 2 and b(k) — \b a g (k)\ (neglecting 
the filter function). Any stochasticity (represented by the 
higher order corrections in P p ) will degrade the correlation, 
yielding different estimates for the bias. Because of this, we 
will consider both the cross- and auto-spectra to check for a 
potential stochastic element of the bias. Q 

As for the density, we can define two estimates for the 

1 Although the effective first-order bias model by DS is deter- 
ministic in Fourier space, it is stochastic in configuration space 
llMatsubaralll999l, iDesiacques fc Shethlfeoifjl) . 
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velocity bias, that we denote be and be c s, with a correlation 
coefficient re. There is no conclusive way to subtract the 
shot noise for the velocity divergence spectrum. Hence 
we used our simulations to gain some insight into this 
issue. We introduced an artificial shot noise in the matter 
spectrum V by randomly drawing a fraction of the particles 
and found that T-shot(fc) — A ■ k 2 asymptotically for large 
fc. Therefore we fitted the amplitude factor A to obtain 
shot-noise corrected power spectra. Note however that 
other terms could be important at smaller wavenumbers; in 
this work we only consider data for which V s hot(k) < 0.1V. 

The DS model describes the biasing of linear density 
peaks that are expected to form haloes of a given mass 
at a specified redshift according to some collapse model 
(which determines the value for <5 C ). However, we analyse 
the density and velocity fields for the actual FOF and SO 
protohaloes. These are the quantities of physical interest for 
st udyin g galaxy clustering in te r ms of halo occupatio n mod - 
els (e.g. ICoorav fc Shethi 120021 ). iLudlow fc Portianil (|201lh 
showed that the vast majority of haloes in our N-body sim- 
ulations can be unambiguously associated with linear peaks 
in the initial conditions when smoothed on the mass scale 
of the halo. For example, Si 70 per cent of all haloes can 
be matched with similar-mass peaks in 5s (x), and Si 90 
per cent for haloes with M Si 10 14 /i _1 Mq. Note, however, 
that the correspondence between the DS peaks and pro- 
tohaloes will not be perfectQ On a given mass scale, some 
peaks with overdensities at the collapse threshold will evolve 
into substructures contained within larger virialized haloes 
(the so-called cloud- in-cloud problem), or to haloes of signif- 
icantly different mass. On top of this, numerical simulations 
have shown that the collapse threshold S c for haloes of a 
given mass and redshift has a broad probabilistic distribu- 
tion r ather than a fixed val ue l|Porciani. Dekel fc Hoffmanl 
|2002| ; iRobertson et al.l l2009h and possibly also depends on 
the form of the adopted smoothing kernel. We consider some 
of these issues in the following section. 



3.3 Barrier heights for top-hat and Gaussian 
niters 

The peak model described in Section [2] requires well- 
defined spectral moments. However, due to its sharp bound- 
ary in real space, the top-hat filter decays very slowly in 
Fourier space so that the integral defining erf is divergent 
in a LCDM model. Since eqs. (fTT]> . (fT2|) and HU) depend 
on 02, DS instead adopted a Gaussian filter and assumed 
S c = Ssc — 1.68. Here, 5 C corresponds to the critical den- 
sity for the collapse of a spherical top-hat perturbation in 
an otherwise unperturbed EdS universe, and this does not 
necessarily apply to peaks in a smoothed Gaussian random 
field. Another issue is that the validity of the simple spher- 
ical collapse model is questionable, at best; the probabil- 
ity of a protohalo or a peak being spherical is null, since 
it would req uire the three eigenva l ues o f the tidal tensor 
being equal. ISheth. Mo fc Tormenl (j200l[ ) showed that the 



2 The measured abundance of haloes in the simulations is of the 
same order of magnitude as (albeit a bit hi gher than) the numbe r 
of peaks with v > 1.68/<ro computed as in lBardeen et al. I dl98d) . 
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Figure 1. Mass-dependence of the linear overdensities measured 
at protohalo centres in the initial conditions of our simulations. 
The shaded regions and contours show the full halo distribution 
after smoothing with a Gaussian filter; connected circles highlight 
the median trend. All values of 8^ have been linearly extrapolated 
to z = 0. Open and filled points correspond to our 150 /i — 1 Mpc 
and 1200 h" 1 Mpc boxes, respectively. For comparison, we also 
show, using boxes, the median trends obtained after smoothing 
with a top-hat kernel. In both cases, the median trends are well 
described by eq. II 241 : the dot-dashed curve shows the SMT result 
(note that this is not a fit to our simulation data) , and the dashed 
curve shows the result for slightly different values of the free pa- 
rameters: 5* = 1.15, a = 0.37 and f) = 0.515. The top panel shows 
the mass dependence of the scatter about the median trend for 
the Gaussian-filtered case, which is well approximated by a simple 



power-law, S = 0.4 <t ( 



6/5 
' 



barrier height in the m ore general ellipsoidal collapse model 
(|Bond fc Mverslll996t ) can be approximated by 



Sec(M, 



5, < 1 + a 



al{R s {M),z c ) 



8'i 



(24) 



where 8, — 5 SC is taken from the spherical collapse model, 
and a = 0.47 and /3 = 0.615 are determined from fits to 
the model results. The presence of the dispersion, ctq, in eq. 
(|24fl results in a mass-dependent barrier height: lower mass 
haloes require, on average, higher overdensities for collapse 
since they must hold themselves together against larger tidal 
forces. It should be emphasized that eq. (|24|l describes the 
mean barrier height; the scat ter about the mean ca n be ap- 
proximated by E = 0.3 o"o (|Robertson et all l2009h . These 
values are valid only for the top-hat filter. 

What is the appropriate barrier height corresponding to 
a Gaussian filter? In Figure Q] we show the linear overdensi- 
ties measured at protohalo centres of mass after smoothing 
with a Gaussian filter on the halo mass scale. The shaded 
regions show the halo data, and connected circles the me- 
dians of the distribution. Open points correspond to results 
from our 150 h~ l Mpc box simulation, and solid points to 
our 1200 /T 1 Mpc box run. Squares show the median 8 C (M) 
for the same sample of haloes, but after smoothing the lin- 
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ear density field with a top-h at filter instead (note t h e goo d 
agreement with the result of ISheth. Mo fc Tormenl (|200 ll ). 
shown as a dot-dashed line in Figure [T}. All values have been 
linearly extrapolated to z — 0. 

It is clear that adopting a Gaussian filter entails lower 
values of 8 C at all massed. Nonetheless, our results are still 
accurately described by eq. (|24[) , albeit with slightly different 
values for the model parameters: <5„ = 1.15, a — 0.37 and 
/3 — 0.515. We show this explicitly in Figure[T]using a dashed 
line. The one-sigma dispersion about the mean trend is well 
described by a simple power-law, E = 0.4 Og 5 , as seen in 
the upper panel of the plot. 



4 RESULTS 

In this section we test the predictions of the DS model 
against the measured density and velocity bias of dark mat- 
ter protohaloes. In order to do so, we calculate the (model- 
predicted) bias parameters b v , b( and b a for each mass bin in 
both of our simulations. Since the values of b v and b^ depend 
on peak height (and hence smoothing scale) we adopt two 
models for the collapse barrier: one assumes the spherical 
collapse model (hereafter SG) and the other the ellipsoidal 
collapse model (hereafter EG). A Gaussian filter is used in 
both cases. Table 1 lists the bias parameters predicted by 
the DS model for these choices of collapse barrier. 

4.1 Density spectra 

In Figure[5]we plot b(k) and b e g(k) extracted from our sim- 
ulations(3 along with the model prediction (eq. (|17|) V Al- 
though the model is not expected to work for very small 
scales, we nonetheless show the results for each box up to 
k ~ 1/Rf (we remind that R{ is the smoothing scale needed 
for the velocity field). Overall, the model expression for the 
initial peak bias is able to describe the simulation results rea- 
sonably well. This can be seen from the solid lines in Figure 
[2] where b v and b^ have been treated as free parameters and 
fitted to both the cross- and the auto spectra. However, it 
is clear that the values for the fit parameters are different 
in the two cases (see Table[2]), implying b(k) 7^ |b c ff(fc)| or, 
equivalently, r(k) ^ 1. In particular, for the range of masses 
in Bin OS, b e ff is negative, i.e. r < 0. The predictions from 
the SG and EG barriers, as apparent from the dashed and 
dotted lines in the right panel in Figure [51 are also negative. 
Hence the bias model matches more closely b e g rather than 
b. This is completely expected, because eq. (|17[) is exact only 
for the cross-correlation between peaks and matter while it 
neglects higher-order corrections for the peak autocorrela- 
tion. However, neither the SG nor the EG barrier provide 
the appropriate values for the coefficients. In particular, the 
EG barrier performs better for low masses (bin OS), while 
the situation is reversed for the higher mass bins. 

3 At any mass scale, M, the smoothing length of a Gaussian filter 
exceeds that of a top-hat filter by a factor of (3-^/ 7I-/2) 1 / 3 Si 1.55. 

4 The values directly measured from the simulations are rescaled 
by the growth factor D(zi n ) to match the theoretical estimates, 
where 8 is linearly extrapolated to z = 0. Note that the actual 
Lagrangian bias is a factor D~ 1 ~ 40 — 50 larger than the values 
reported in the figure. 



Figure [2] suggests that the stochasticity is more of a 
problem for haloes with M < M,, and on smaller scales. 
Since high-mass haloes are h ighly correlated with de nsity 
peaks in the initial conditions (|Ludlow fe Porcianill201ll ') and 
peaks follow eq. (|14[) . the relationship between <5h(k) and 
5(k) is likely to be more deterministic for massive haloes. 
Note, however, that the estimate of r is affected by the shot- 
noise correction, which is large in our samples. Therefore we 
cannot draw definitive conclusions regarding stochasticity. 

The performance of the model is characterized in an- 
other way in Figure [3l where we compare directly the mass 
dependence of the bias parameters b v {M) and b^(M) pre- 
dicted by the SG and EG models to the best-fit values ob- 
tained using the parameterization of the cross-spectrum in 
eq. (|17p . We do not show the corresponding results from 
the auto-spectrum because of the uncertain shot-noise sub- 
traction. Figure [3] reveals that overall the b v values derived 
from the SG model are closer to the fitted ones. In contrast, 
there is almost no difference between the two models for b^. 
For both parameters the mass dependence follows a different 
trend when we compare the fits to the models. In particular, 
the prediction for b^ and M < M» becomes more and more 
inaccurate with decreasing halo mass (up to a factor of ~ 3 
for M ~ M«/15). On the other hand, in the mass range for 
which the model was developed (M > M„), b^ is overpre- 
dicted by a factor of ~ 1.5, independent of halo mass. Figure 
shows that the disagreement is even worse when b v and 6^ 
are obtained from P^. All of this demonstrates that eqs. ()11|) 
and (|12|) cannot accurately describe the mass dependence of 
the Lagrangian density bias, especially for M < M 4 . We will 
discuss possible reasons for this in Section O but first turn 
our attention to the velocity spectra. 

4.2 Velocity spectra 

Figure [4] shows be(k) and be c s(k) for bins 1S-4S and 1L-4L; 
the theoretical model, b ve i(k), is also plotted. The data have 
been "de-smoothed", i.e. they have been divided by a filter 
function with smoothing scale Rf, the same used to smooth 
the velocity field in the first place. It is apparent that on 
large scales the two estimates of the bias coincide, i.e. re ~ 1, 
indicating a strong correlation between the fields. Overall, 
we can conclude that the velocity bias is deterministic to a 
good approximation. 

There is an excellent agreement (better than 10 per 
cent) between the simulation results and the model predic- 
tions for all mass bins at scales k < 0.1 /iMpc -1 for the 
large box and k < 0.3 h Mpc - for the small one. Therefore, 
the peak model provides a faithful description of the velocity 
bias. 



5 SUMMARY 

We have investigated the Lagrangian bias of dark-matter 
haloes by testing the theoretical model proposed by DS 
against N-body simulations of structure formation. The 
model assumes that haloes form from density peaks and pre- 
dicts a scale-dependent bias for both the density and velocity 
fields. Our main results can be summarized as follows. 

• When averaged over a spherical Lagrangian volume 
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Table 1. Mass range, number of haloes, average mass and bias parameters (both weighted by halo counts) for the nine different mass 
bins used in our study. The first five belong to the small-box simulation; the last four to the large-box simulation. The bias parameters 
are computed according to the two models for the peak height described in the text: the spherical and ellipsoidal collapse models. All 
values have been obtained adopting a Gaussian filter. 
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Figure 2. Left: Stars and open circles plot, respectively, b(k) and b c g(k) for the four separate mass bins in our large simulation box. 
Dotted lines show the predictions obtained for the spherical collapse model, and dashed lines those of the ellipsoidal collapse model 
adopting a Gaussian filter. Solid lines are the best-fits to the data points. Right: Stars and open circles correspond to b(k) and \b c g (k)\, 
respectively. As in the left panels, dotted and dashed line show to the SG and EG models, and solid lines are the best-fits to the data. 
Thick lines indicate negative values of the bias. In all cases, errors are propagated assuming that the uncertainty in the power-spectrum 
is c(fc) = 2/N(k) ■ P(k), where N(k) is the number of modes in each k bin. 
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Table 2. Best-fit values for the density bias parameters from the cross-spectra and the auto-spectra. Uncertainties are always at the few 
per cent level. Haloes are identified cither with the FOF or the SO algorithms. 
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Figure 3. Density bias parameters, b v (left) and 6^ (right), plotted as a function of halo mass. Orange circles with error bars correspond 
to the best-fit values obtained from the cross-spectra of haloes identified with the FOF algorithm. The corresponding values for SO haloes 
are shown with purple stars. The solid and dashed linos show, respectively, the predictions of the DS model for the spherical collapse 
model, and for the ellipsoidal model with a Gaussian filter. The vertical dashed black lines mark Af*. 
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Figure 4. Stars and circles show bg(k) and bg a g(k), respectively, for haloes in our large (left) and small (right) simulation boxes. The 
solid line is the DS model prediction (eq. H15H with R s given in Table [TJ, which is independent of peak height. Errors are propagated as 
in Fig. [2] Note that we only show bg(k) on scales over which the shot-noise is sub-dominant, as stated in Section 13.21 



containing the appropriate mass, the linear density contrast 
measured at protohalo centres depends sensitively on the 
choice of smoothing kernel. For a Gaussian kernel, for exam- 
ple, the resulting density contrasts are systematically lower 
than those computed using a top-hat filter. This is because, 
at fixed mass, the smoothing length of a Gaussian kernel 
exceeds that of a top-hat filter by a factor of about 1.55 re- 
sulting in systematically lower density estimates. Nonethe- 
less, the median barrier height computed with a Gaussian 
kernel can be accurately p arameterized by the same fit ting 
formula first advocated bv lSheth. Mo fc Tormenl ^OOll ) for 



the case of a top-hat filter, albeit with different values for 
the numerical parameters. We use this result to approximate 
the collapse threshold for dark matter halo formation when 
adopting a Gaussian filter. 

• The functional forms for the density and velocity bias 
relations derived by DS - our eqs. (|17|) and (|18[) - accurately 
describe the results obtained from our simulations, provided 
the parameters of the model are allowed to vary with respect 
to the model-predicted values. In both cases, the Lagrangian 
bias is characterized by a constant term that dominates on 
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very large scales, and a scale-dependent term proportional 
to k 2 . 

• Quantitatively, the velocity bias predicted by the DS 
model is able to reproduce the measured protohalo velocity 
bias in our simulations to better than 10 per cent, provided 
we limit ourselves to quasi-linear scales (fc ^ 0.3 /tMpc -1 ). 
The predicted density bias (eqs. and (H2J), on the other 
hand, does not match the Lagrangian density bias extracted 
from the simulations. This is likely due to the more com- 
plex nature of the density bias, which depends additionally 
on peak height and on the exact definition of a protohalo. 
These results are independent of whether one adopts a bar- 
rier height consistent with either the spherical or ellipsoidal 
collapse model. 

• We have measured the mass dependence of the density 
bias, b v and b(, by fitting the density power-spectra obtained 
for haloes in several different mass bins in each of our two 
simulations. The most massive haloes identified in our sim- 
ulations are the most strongly biased, and have characteris- 
tic overdensities that correlate with the underlying matter 
density. In contrast, the distribution for low mass haloes 
(Mh < 10 12 /i _1 Mq) exhibits an anti-correlation. A similar 
trend is also predicted by the DS model, although with no- 
ticeable differences (see Figure [3]). We emphasize that the 
best-fit values for b v and b^ obtained in our analysis apply 
for a LCDM cosmogony (with our adopted cosmological pa- 
rameters) and only over the limited mass range probed by 
our simulations. Future work should consider how the La- 
grangian density bias depends on the underlying cosmology. 

• Our comparison of the auto- and cross-spectra in Sec- 
tion|3]suggests the presence of a stochasticity in the Fourier- 
space density bias of dark matter protohaloes, but none 
for the velocity bias. This is in disagreement with the DS 
model for the bias at leading order, and seems to cor- 
roborate the need for higher-order terms in the expression 
for the auto-spectrum of the protohaloes, as predicted in 
iDesiacques et ail (|2010T ). However, we cannot draw firm con- 
clusions regarding stochasticity in the bias estimates due to 
uncertainties in the shot-noise subtraction. 

We have tested the Lagrangian bias model of DS against 
a pair of high resolution N-body simulations of structure for- 
mation and found that it is able to accurately reproduce 
the velocity bias, but not that of the density. Our anal- 
ysis focused on protohaloes, the high redshift progenitors 
of 2 = dark-matter haloes, whereas the model describes 
the biasing of density peaks. One possible explanation for 
the differences in the model predictions and simulation re- 
sults comes form the differences between the expressions for 
the bias parameters: 6„, fe^ and b a . The latter of the three 
is solely determined by the linear matter power spectrum 
and the smoothing scale corresponding to a given mass, M. 
The density bias parameters, however, depend additionally 
on explicit properties of the peaks, such as their height, v, 
mean curvature, u, and on an assumed collapse threshold for 
their identification as haloes of a given mass. Figure [T] shows 
that there is a large halo-to-halo variation in v at any given 
mass scale; characterizing the collapse threshold as a single 
mass-dependent value may, therefore, be too simplistic. This 
added complexity introduces a significant margin for error 
in the model's estimates of the density bias. It remains to 
be seen whether more realistic models of the collapse bar- 



rier - such as those that account for the statistical scatter 
in the linear over-densities of protohaloes at a given mass - 
will improve the model's predictive power. 

Another possibility for the discrepancies stems from the 
fact that we are identifying linear density peaks with pro- 
tohaloes in our simulation initial conditions. Although the 
majority of our dark matter protohaloes form in the vicin- 
ity of linear density peaks of the same characteristic mass, 
the fate of all peaks with the same mass and overdensity is 
unclear. Uncertainties associated with the identification of 
protohaloes within the linear density field may have adverse 
effects on the predictive power of the density bias model. A 
more detailed understanding of how protohaloes map onto 
linear density peaks, and vice versa, will, no doubt, provide 
valuable insight into the mechanisms behind halo biasing. 
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